function [S1minus,S1plus,S2minus,S2plus,S3minus,S3plus,dt1,dt2] = BiEcllipticNonPlanar(R1,R3,delta_i,r2_ratio,mu);

%% COE of initial state

COEi = [norm(R1),0,delta_i,0,0,0]';
%[] Initial state

S1minus = Coe2State(COEi,mu);
%[] initial orbit's state before delV

%% COE after first burn

r1 = norm(S1minus(1:3));
%[] Range of the initial orbit

r2 = r1 * r2_ratio;
%[] Range of the transfer orbit 1

sigma = r1 / r2;
%[]Ratio of orbital ranges.

emin1 = (1 - sigma) / (1 + sigma);
%[]Eccentricity.

COE_T1f = [(COEi(1)+r2_ratio*COEi(1))/2,emin1,COEi(3),0,0,0]';
%[] COE of the spacecraft after first burn

S1plus = Coe2State(COE_T1f,mu);
dt1 = pi*sqrt(COE_T1f(1)^3/mu);
%% COE before second burn

COE_T2o = [(COEi(1)+r2_ratio*COEi(1))/2,emin1,COEi(3),0,0,180]';
%[] COE of the spacecraft after first burn

S2minus = Coe2State(COE_T2o,mu);


%% COE after second burn
COEf = [norm(R3),0,0,0,0,0]';
r3 = COEf(1);
sigma2 = r2/r3;
emin2 = (1 - sigma2) / (1 + sigma2);
COE_T2f = [(COEf(1)+r2_ratio*COEi(1))/2,-emin2,COEf(3),0,0,180]';
S2plus = Coe2State(COE_T2f,mu);
dt2 = pi*sqrt(COE_T2f(1)^3/mu);
%% COE before third burn
COE_T3o = [(COEf(1)+r2_ratio*COEi(1))/2,-emin2,COEf(3),0,0,0]';

S3minus = Coe2State(COE_T3o,mu);

%% Coe after third burn

COE_T3f = COEf;

S3plus = Coe2State(COE_T3f,mu);

end